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The advisability of constructing a CAD system for glass-melting furnaces is demonstrated. The main princi¬ 
ples of the mathematical support for the CAD technological subsystem are formulated. The principles of con¬ 
structing a mathematical model of the thermal performance of a furnace are described. The main equations 
providing for a combined solution of the exterior and the interior problems of heat and mass exchange taking 
into account the specifics of the physiochemical aspects of the glass-melting process are given. 


The efficiency of investment projects in the glass indus¬ 
try to a great extent depends on the technical and economical 
parameters of the particular glass-melting furnace. The high 
capital cost of a modern furnace, which is primarily deter¬ 
mined by the cost of refractory materials, can be justified in 
the case of high efficiency, low service cost, and a protracted 
furnace campaign [1]. The cost of constructing a furnace in 
Russia approaches these costs in industrial countries. Conse¬ 
quently, an investor has the right to expect a corresponding 
technical efficiency from the furnace. In this context the 
quality of the furnace design becomes critical. 

It should be admitted that designing a furnace is tradi¬ 
tionally reduced to the determination of geometrical sizes of 
particular furnace components. The design calculations at 
best include the determination of the thermal balance items 
for the melting and working zones of the furnace, the volume 
of the checkerwork regeneration, and other local calculations 
(parameters of the smoke funnel, burners, etc.). This scope of 
calculations can suffice only for a draft design of the furnace, 
which traces its layout, the positions of its fundamental sup¬ 
ports, and the loads on them. The detail design of the furnace 
implies a deeper and more comprehensive account of all 
physicochemical aspects of the glass-melting process. 

Regarding a glass-melting furnace as a continuous en¬ 
ergy-technological melting aggregate, one can isolate two in¬ 
terrelated sets of processes that determine the design and ef¬ 
ficiency of the furnace. The first set consists of generation of 
thermal energy in the working space and in the tank (addi¬ 
tional electric heating); the second one implies the transfer of 
heat to the process zone. The thermal stress value and its dis- 
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tribution over the flame space of the furnace are determined 
by the fuel combustion intensity and gas dynamics of the 
combustion products flow, i.e., the set of the radiation and 
aerodynamic characteristics of the flame. Other terms being 
equal, the length of the flame, its velocity, and its position 
with respect to the heated surface have a determining effect 
on temperature distribution over the glass melt surface and 
the enclosing refractory brickwork. The above parameters fi¬ 
nally detennine the efficiency of exterior heat exchange and 
the integral heat absorption in the tank [2]. 

The distribution of the heat absorption of the melt in the 
volume of the process zone to a great extent depends on the 
thermophysical characteristics of the melt and the hydrody¬ 
namics of the tank. The intensity and nature of convection 
glass flows depend not only on the temperature field of the 
glass melt surface, but also on the design of the melting tank 
and the location of the bubbling system and barrier devices. 
When additional electric heating is used, the location and 
power of additional thermal power sources have a fundamen¬ 
tal importance for the heat absorption and hydrodynamics of 
the melting tank. 

Thus, at the stage of a detail design, which is the second 
design phase, it is necessary to obtain full quantitative infor¬ 
mation on the temperature fields of all furnace elements par¬ 
ticipating in heat exchange. This includes the temperature 
distribution on the roof and on the lateral walls of the upper 
furnace structure, the melt surface, the walls and the bottom 
of the melting tank, and the gaseous space of the furnace and 
the melt. This information assists in making a justified deci¬ 
sion on the optimal characteristics of the flame, the choice of 
refractory materials for brickwork elements that undergo the 
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most substantial wear, identifying stagnation zones in the 
melting tank, determination of the final layout of the barrier 
devices, and, if necessary, the input site and power of addi¬ 
tional electric heating. Such information, as a rule, promotes 
a reconsideration of the initial geometric layout of the fur¬ 
nace. 

Based on the above, it can be stated that a correct design 
of a glass-melting furnace can be implemented only applying 
an automated design system, of which CAD is the optimum 
form. The development of CAD for a glass-melting furnace 
is a most complicated problem. The first stage of solving this 
problem involves the development of a number of local sub¬ 
systems to provide solutions for geometrical and technologi¬ 
cal design problems. The most complicated subsystem is the 
technological design subsystem, which cannot be created 
without corresponding technical, mathematical, software, 
and informational support. 

Analysis indicates that the technical support of a CAD 
subsystem includes an up-to-date PC with an extended pe¬ 
ripheral infrastructure. The majority of the design proce¬ 
dures, including the construction of functional models, geo¬ 
metrical design problems, and numerical methods for solving 
algebraic and differential equations are implemented using 
standard software products. Published results of theoretical 
and experimental studies of glass melting and the 
thermophysical properties of batches, melts and refractories 
provide a database that can be used as general informative 
support. Thus, the main difficulty in developing a technolog¬ 
ical design subsystem involves the development of its mathe¬ 
matical and special informational support based on a mathe¬ 
matical model of the thermal performance of a glass-melting 
furnace. Let us discuss the main concepts of this model. 

Three groups of equations were used to formalize heat 
and mass exchange processes in a glass-melting furnace: 

- exterior problem equations describing the regularities 
of heat release and mass transfer in fuel combustion, radiant 
and convective heat exchange between hot and burning com¬ 
bustion products and surfaces of the tank and enclosing 
brickwork of the working space; 

- interior problem equations determining the propaga¬ 
tion of heat received by the tank surface over the tank vo¬ 
lume, taking into account the kinetics of physicochemical 
transformations and heat losses via the tank brickwork; 

- tank hydrodynamics equations describing the effect of 
the parameters of the exterior and interior problems on the 
three-dimensional convection of the glass melt. 

In developing the numerical model of the furnace, the 
main system theory requirement (rule) was taken into ac¬ 
count, namely, that all processes taking place in a object con¬ 
stitute an integrated system and are interrelated and mutually 
determined [3], 

The resolvent zonal method [3, 4] was used for solving 
the exterior problems. This method implies splitting the 
melting zone of the glass-melting furnace into volume and 
surface zones with preset geometrical sizes. The gas spaces 


is partitioned into volume zone, whereas the enclosing brick¬ 
work and the glass melt mirror are split into surface zones. 

The heat exchange energy characteristics for each partic¬ 
ular zone are determined by solving a system of nonlinear al¬ 
gebraic equations of thermal balance and heat transfer writ¬ 
ten in the following form: 

n + m n+ m 

X a ji T i + YjSji T i + Qi = 0, 1 = 1,..., m + n, (1) 

j =1 7=1 

where m and n are the quantities of the volume and surface 
zones, respectively; a- is the coefficient of radiation exchange 
between zones j and i; 7j is the temperature of zone i, K; g.. is 
the coefficient of mass exchange and convective exchange 
between zones j and i; Q ; is a free summand of the equation 
taking into account the exterior and interior heat sources and 
outlets. 

The radiation exchange coefficients characterizing the 
quantity of radiation energy transferred from zone j to zone i 
at T. = 1 are found from the equations: 

for volume zones 

a ji 4 Vj ^0 kjfji ; 

for surface zones 

a ji — f / C-0 > 

where V, and F. are the volume and the surface area of the re¬ 
spective zones; c 0 is the Stefan - Boltzmann constant, 
c 0 = 56.687 x 10- 12 kW/(m 2 • K); k j is the reduction coeffi¬ 
cient of zone j, m“f. is the resolving generalized angle fac¬ 
tor characterizing the share of energy absorbed by zone i 
from zone j taking into account that part of the energy of the 
radiation ray is absorbed by gas; e:■ is the degree of blackness 
of the surface zone j. 

The coefficients / /7 are found by solving a system of li¬ 
near algebraic equations using the algorithms proposed by 
V. G. Lisienko and Yu. A. Zhuravlev [5,6]: 

for volume zones 

n 

fji='Yji + 'L r kVjkfki ; ( 2 ) 

k = \ 

for surface zones 

n 

fji = ' Vji + YaWjkfki » Uj =1,2,..., m + n, (3) 

*=i 

where \|/ 7 is the generalized angle coefficient of the radiation 
from zone j to zone i; \\i- k is the same for radiation transferred 
from a surface zone, for instance, k\ r k = 1 - s k is the reflect¬ 
ing capacity of the zone k. 

In Eqs. (2), (3) the first summand in the right-hand side 
characterizes the direct heat transfer from zone j to zone i. 
The second summand describes the heat transfer between 
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these zones taking into account the reflection of radiation 
from other surface zones k that are “visible” for zones j and i. 

The coefficients \|/.. are calculated in accordance with 
the statistical test method using the algorithm proposed by 
A. E. Klekl’ and S. D. Dreizin-Dubchenko [7, 8]. Each fixed 
element i in this set corresponds to a system of m + n linear 
equations. The z'th column of the matrix \|/ is found on the 
right-hand side of this system. The unknown quantities are 
the elements of the zth column of the matrix f. Thus, to calcu¬ 
late the resolving angle coefficients with the constraint 
i = 1, ... m + n, it is necessary to solve m + n systems of lin¬ 
ear equations that differ only in their right-hand sides. This is 
implemented using a standard subprogram in Compaq Visual 
Fortran 6.6.B Professional. 

The summand Q t depends on the specific statement of 
the problem and includes heat release in fuel combustion, the 
heat consumed to compensate the thermal effect of chemical 
reactions, the heat used for heating of the batch, cullet, and 
the glass melt, the heat losses via the walls and bottom of the 
furnace, etc. 

The system of equations (1) has many parameters whose 
values are determined by iteration. This primarily concerns 
the radiation characteristics of combustion products, which 
depend on the temperature and composition of the gas. 
Therefore, a combined approach is used for modeling the 
thermal work of a glass-melting furnace. Initially the temper¬ 
ature of the combustion products and the degree of burn-out 
of fuel in particular zones of the working space are set using 
the typical conditions of the operation. According to the pre¬ 
set values of the degree of combustion, the composition of 
the gaseous phase is determined and the required characteris¬ 
tics are calculated using the standard method [9], The proce¬ 
dure described is repeated in the context of the zonal method 
until the convergence of the results. 

Nevertheless, the calculation results in such an approxi¬ 
mation cannot be regarded as final, since the zonal method 
does not take into account the regularities of the gas flow and 
fuel combustion under this flow. To eliminate this drawback, 
the zonal-node method is used [10]. The temperature fields 
and the combustion product fields determined by the zonal 
method are approximated by continuous functions using the 
Table Curve 3D software package. The working space of the 
furnace is covered by a nonuniform finite-difference grid in 
such a way that the nodes of the volume (gaseous) zones co¬ 
incide with the grid nodes. The flow and heat and mass trans¬ 
fer equations for the working space of the furnace are solved 
using the SIMPLER method developed by S. Patankar [11] 
using the approximation results as initial approximation. 

The most general form of gas mechanics and heat and 
mass transfer equations taking into account the combustion 
of fuel is presented in [12, 13]. The mathematical model of 
the glass-melting furnace took into account the specifics of 
the thermal performance of the particular furnace related to 
the cross-link effects, i.e., the mutual influence of specific 
phenomena on each other. Thus, it was taken into account 


that the mass transfer caused by barodiffusion with respect to 
the working space of the glass-melting furnace can be ne¬ 
glected due to the low pressure gradient as compared to the 
pressure itself. Within the temperature interval up to 2500 K 
the values of the thermal diffusion coefficient are within the 
interval of 10 “ 3 - 10 5 . Therefore, although the temperature 
gradient at the combustion front may be substantial, the 
mechanism of thermal diffusion in general can be neglected 
as well. Thus, the equations of the “fuel” and “oxidizer” ba¬ 
lance in the mathematical model of the glass-melting furnace 
have the following form: 

(Qyvi A 

pi ^ + V grad m f I = div (pD ef grad m f ) + R f ; (4) 

(dm ^ 

P i + V S1 " ad “ 0X J = ^ (PDe£ ‘ gI " ad m ° X ' >+R ° x ’ (5) 

where p is the density of the mixture of gas, air, and combus¬ 
tion products, kg/m 3 ; m { and m are the mass contents of the 
“fuel” and the “oxidizer”; D eU is the effective diffusion coef¬ 
ficient of the zth component in the mixture, m 2 /sec; R { and 
R ox are the flow rates of the “fuel” and the “oxidizer”, 
kg/(sec • m 3 ); v is the mass center velocity vector with com¬ 
ponents Vj = v x , v 2 = v , v 3 = v z , m/sec; z i is the generalized 
coordinate with the components Xj =x, x 2 = y, and x 3 = z 
(i.e., i ',7 = 1,2,3), m. 

For a simple chemical reaction where two components 
(the fuel and the oxidizer) react in a certain ratio [1 kg 
“fuel”+ z'kg “oxidizer” —> (l+z')kg reaction products], 
their flow rates are found from the expressions 

Rf—Kp 2 m { m ox exp 



R ox = iR { , (6), 

where K, E are empirical coefficients determined by the kind 
of fuel and by the type of its combustion. 

The closure of the gas mass balance for the working 
space is implemented taking into account the equation of the 
total gas (mixture) continuity: 

^ + div (pv) = 0. (7) 

ox 


Using the mass center velocities v, the equations of com¬ 
pressible gas flow with variable viscosity have the form of 
the traditional Navier equations: 


dv 

p— = -gradp + 
dr 


Div 


bef 


& 


8Z: 


2 

3 


p ef (divv)5 



( 8 ) 


where p^. is the partial density of the Mi component of the 
mixture of gases, kg/m 3 ; x is the time, sec; p is the mixture 
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pressure, Pa; p ef is the effective dynamic viscosity of the 
mixture, Pa • sec; 5 is the identity second-rank tensor, whose 
elements are Kroneker symbols 8 (/ ; g is the acceleration vec¬ 
tor of the gravity per mass unit of the component k\ Div is the 
stress tensor divergence; div is the vector divergence; dv/dx 
is the total or substantial velocity derivative (acceleration) of 
the mixture determined by the relation 


The equation of the transfer of turbulence energy dissipa¬ 
tion velocity is 


dz _ 

-H v grade = 

St 


div 


7 > 


V t 

V H-— 

[l Pl «v 

grade 


+ C, —-C 


E t E t 


( 12 ) 


dv 

dx 


dv 

+ (V V) v, 
St 


where V is the vector differential operator “nabla” [3]. 

The gas flow in the working space of the glass-melting 
furnace has a clearly expressed turbulent nature, hence, the 
effective coefficients p ef = p + p t and £> ef depend on the de¬ 
gree of the flow turbulence. There are several engineering 
theories of turbulence known. The most convenient in the 
practical aspect is the Kolmogorov - Prandtl model E t — s 
[3, 14], which has been used in our study. According to this 
model, the turbulent component of the effective dynamic vis¬ 
cosity |_i t is proportional to the square of kinetic turbidence 
energy E { and inversely proportional to the velocity of the 
turbulence energy dissipation s. Then the kinematic turbu¬ 
lence viscosity v t is determined by the expression 


v _ Pt _ C » E ? 

v t - 

P s 


(9) 


where C is the empirical coefficient; C p = 0.09; E { is the ki¬ 
netic turbulence energy per weight unit: 

£, = 0.5 [((v;) 2 ) + ((v;) 2 ) + <(v;) 2 )]; 


here the angular brackets denote the time averaging and the 
primes indicate the pulsation constituents (fluctuations) of 
the velocity components. 

Specific balance equations (transfer equation) are imple¬ 
mented for the variables E t and s. 

The equation of kinetic turbulence energy transfer is 


8E t 

dx 


■ v grad E t = div 


Pr, 


gradfs] 


t 


+ G -8, 


( 10 ) 


where C j, C 2 , Pr t v and Pr E are empirical coefficients taken 
from the intervals Pr t v = 0.9 - 1.0, Pr E = 1.22 - 1.3, C] = 
1.44- 1.55, and C 2 = 1.92-2.0. 

The gas flow energy equation (the Fourier-Kirchhoff 
heat transfer equation) takes into account all the main and ac¬ 
cessory (cross-link) effects of its enthalpy variation [3, 12]. 
The summands on the right-hand side of the equation charac¬ 
terize heat transfer by radiation, heat conduction, and diffu¬ 
sion heat conduction, as well as heat transfer due to the pres¬ 
sure gradient; heat sources (outlets) caused by chemical reac¬ 
tions; work of pressure forces; heat sources cause by flow 
energy dissipation, i.e., by the conversion of the internal fric¬ 
tion work into heat; heat transfer due to the diffusion mass 
flow and heats sources caused by the work of exterior forces 
in the diffusion mass flow process. 

Some enthalpy variation sources, as applied to the condi¬ 
tions of the glass-melting furnace, can be neglected in view 
of their small values. It has been noted that the pressure gra¬ 
dient is very small; consequently, the thermal effect of pres¬ 
sure can be neglected. The work of pressure forces and aero¬ 
dynamic heating is not significant as well. Heat transfer un¬ 
der the effect of gravity and the diffusion thermal effect are 
extremely insignificant. Consideration of the above circum¬ 
stances allows for a perceptible simplification of the equa¬ 
tion of the energy of the medium. At the same time, addi¬ 
tional aspects should be taken into account. This primarily 
concerns the thermal conductivity of the glass melt, which 
should be replaced by its effective value A, ef equal to the 
sum of its conductive and radiation components. Next, one 
should take into account the value E t . Furthermore, the Fou¬ 
rier - Kirchhoff equation should be transformed for turbulent 
flow conditions. 

Combining the energy equation with the equation of the 
gas mixture continuity (7), we can write the equation of the 
thermal balance of gases suitable for the analysis of heat 
transfer in the working space of the glass-melting furnace: 


where v = p/p is the molecular kinetic viscosity; G is the tur¬ 
bulence generation. 

For the calculation of G we use the expression 


G = v t 




!*lT +2 &] 2 + 

dx J { dy J 




dV 

dy , 



( 11 ) 


—(pT) + div(pvr) 

OX 


= div (A, rf grad T+q R )~Y. H j R j 


j =f 


div 


Fef 


g ra d — + r* grad£ t +p jgradwi • 


j =f 




(13) 


j =f 


where T is the gas temperature, K; c is the specific heat ca¬ 
pacity of the gas at constant pressure, J/(kg • K); q R is the ra- 
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diation flux density, W/m 2 ; A. f is the thermal conductivity of 
gases, W/(m • K); Q* is the heat transfer of the yth compo¬ 
nent of the gas mixture (the heat transferred by one mole of 
this component as a consequence of diffusion at grad T= 0), 
J/kg; Hj is the thermal effect from the reaction of the con¬ 
sumption of the yth component, J/kg; J* is the diffusion flux 

density, kg/(m 2 • sec); = v + v t /Pr t v is the transfer (diffu¬ 
sion) coefficient of turbulence energy, m 2 /sec. 

Note that Eq. (13) is integrodifferential, since the vector 
of the radiation flux density q R has the following compo¬ 
nents: 




| IrijdQ ., 

471 


where I is the radiation flux intensity; n i is the direction co¬ 
sine of the thermal ray with respect to the zth axis of the coor¬ 
dinate system; dfl is the elementary body angle. 

The system of equations (4)-(13) is solved under the 
boundary conditions of the 3rd kinds on the walls and the 
roof of the furnaces and under the boundary conditions of the 
4th kinds (conjugacy conditions) on the surface of the melt¬ 
ing tank. The dynamics of fuel combustion is calculated in 
accordance with Eqs. (4) - (7) in the context of the Schwab - 
Zel’dovich approximation [11]. 

The results of cyclic solving of zonal method equations 
(1) - (3) and gas dynamics and heat and mass transfer equa¬ 
tions (4)-(13) are consistent fields of gas temperatures in 
the working space of the furnace and on the enclosing sur¬ 
faces, gas composition, velocities, and pressures. This makes 
it possible to relate the zonal method to the regularities of gas 
flow, fuel combustion, and heat exchange in the working 
space volume, approach the problem of optimizing the de¬ 
sign and regime parameters of the fuel-burning devices 
(burner nozzle diameters, layout of the burners, selection of 
the angle of the burner with respect to the tank surface, etc.), 
and finally optimizing the geometrical parameters of the fur¬ 
nace. 

It should be noted that the above cited exterior problem 
equations form an interrelated mutually determined chain. In 
principle they characterize the energy aspects of the glass¬ 
melting process. In this sense the concluding equations of the 
mathematical model, i.e., the equations of the tank hydrody¬ 
namics, have a certain independence. Their solution requires 
only information on heat flux densities on the glass melt sur¬ 
face and temperatures on this surface from the exterior prob¬ 
lem. Although these equations are equations of convective 
heat exchange between the melt and the enclosing surfaces, 
i.e., they describe thermal energy transfer regularities, the 
purpose of this block of the mathematical model is purely 
technological. It is intended to identify the thermal condi¬ 
tions of the furnace providing for an optimum homogeneity 
of the glass melt. 


The tank hydrodynamics equations include: 
glass melt flow equations 


gm 


d/grrp 

dx 


+ ( v gm V ) V gm 


grad p gm + Div (^ grad v gm ) + p Em g; 


gm 


gmi 


(14) 


the equation of glass melt continuity 

3p g m 


dx 


- + div (p gm v gm ) = 0; 


(15) 


the equation of glass melt energy combined with the con¬ 
tinuity equation 


(Pgin^gm ) + div(Pgm Vgm^gm ) 
div (Xer.gm grad T gm ), 


(16) 


where p gm is the glass melt density, kg/m 3 ; v gm is the glass 

melt flow velocity, m/sec; p gm is the glass melt pressure, Pa; 
Pgm is the glass melt dynamic viscosity, Pa • sec; c is the 
specific heat capacity of the glass melt, J/(kg • K); T gm is the 
glass melt temperature, K. 

The equations of the discrete analog of system (14) - (16) 
in the three-dimensional version are solved by the SIMPLER 
method; two-dimensional approximation is implemented in 
the “vorticity - current function” variables [3, 14]. Since the 
components of the glass melt velocity on a free surface (open 
channel) are required in any case, the “creep flow” problem 
of glass melt in the tank is primarily solved. Note that the re¬ 
liability of solving the tank hydrodynamics equations to a 
great extent depends on the reliability of data on the thermo¬ 
physical properties of the glass melt and the dynamics of 
physiochemical processes occurring in its formation and 
melting. 

To conclude, it should be noted that the mathematical 
model of the thermal performance of a glass-melting furnace 
considered here enables one to study the regularities of the 
exterior and interior heat and mass transfer. The analysis of 
these regularities will make it possible to solve a number of 
problems in the optimization of different functional relation¬ 
ships and thus to create a special informational support for 
the subsystem of the technological design of glass-melting 
furnaces. 
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